run_therMizer

Purpose
Load libraries
Set up multispecies model
Load data needed for FishMIP scenarios
Write parameters and functions for therMizer extension
Run the simulations

Purpose

The purpose of this script is to run the FishMIP simulations using therMizer. These runs are for the regional model of Hawaii’s longline fishery. They include temperature and plankton forcing from the CMIP6 models. Many thanks to Gustav Delius for several tips on improving earlier versions of the code used here.

Load libraries

Mizer version 2.0.3 or later is needed to use the extension.

if (!require("mizer", character.only = TRUE)) {
  install.packages("mizer")
} else if (packageVersion("mizer") < "2.0.3") {
  install.packages("mizer")
}

library("mizer")

Set up multispecies model

For the FishMIP simulations, we’ll be using the species parameters in Woodworth-Jefcoats et al. 2019. We’ll also follow their implemetation of fishing mortality, so we’ll need to create a selectivity function which we’ll call knife_edge_phased. We’re doing this before creating the params object so that the fishing selectivity parameter is recognized.

In order to use the therMizer extension, two additional species parameters are needed: temp_min and temp_max. These represent the thermal tolerance limits for each species. You could find these in the literature in physiological or tagging studies. Alternatively, they can be assumed based on species’ vertical and geographic ranges, that of their prey species, or some other information. Note that these are only needed for the therMizer extention functions that relate to temperature, not those that relate to plankton.

# Create the selectivity function:

#' Size based knife-edge selectivity function that is phased in 
#'
#' A knife-edge selectivity function where only sizes greater or equal to
#' \code{knife_edge_size1} are selected.  Selectivity is then phased in through
#' size \code{knife_edge_size2}.  This simulates the reduced effectiveness of 
#' longline gear for smaller fish
#'
#' @param w The size of the individual.
#' @param knife_edge_size1 The size at which the knife-edge is initiated.
#' @param knife_edge_size2 The size at which selectivity is fully realized.
#' @export
knife_edge_phased <- function(w, knife_edge_size1, knife_edge_size2, ...) { 
    sel <- rep(0, length(w)) 
    
    # Phased in linearly
    F0 <- which(w < knife_edge_size1) # to find one size smaller than that fished, for the 0 value 
    F1 <- which(w < knife_edge_size2) # to find end of escalation size range 
    lo_sel <- max(F0):max(F1) 
    sel[lo_sel] <- seq(0, 1, length = length(lo_sel)) # linear increase from 0 to F 
    sel[w >= knife_edge_size2] <- 1 
    return(sel) 
    
} 

# Load species parameters
species_params <- read.csv('NPac_species_params.csv')

# Load interaction matrix
inter <- read.csv('inter_NPAC.csv', row.names = 1)
inter <- as(inter, "matrix")

# Create the params object
params <- newMultispeciesParams(species_params, interaction = inter, min_w_pp = 1e-14, no_w = 100, kappa = 1e12, w_pp_cutoff = 455400*1.1)

Load data needed for FishMIP scenarios

This round of FishMIP simulations includes three forcing inputs: fishing mortality, temperature, plankton. We’ll load them in that order.

Fishing

There are two fishing scenarios:

histsoc: Historical effort until 2015, then fixed at 2015 levels
nat: No fishing

The time series of fishing mortality, “FishingEffort.dat” was created with “CreateF.Rmd”. We’re also going to create another time series with F = 0 so that we can use the same years for all our scenarios. This makes incorporating the temperature functions easier later on because all the time values will agree.

# Load fishing scenario
effort_Fhistsoc <- read.table("FishingEffort.dat", sep = " ")
effort_Fhistsoc <- effort_Fhistsoc[,2] # Column 1 = year, which we won't need
effort_Fhistsoc <- as(effort_Fhistsoc, "matrix")

# Build fishing effort arrays
# Since we're following the approach of Woodworth-Jefcoats et al. 2019, we'll begin the model with 600 years of constant input for spinup.  This means the first year will be 1350 (1950 - 600).
times <- seq(1350, 2100, by = 1)
gear_names <- c("Longline")
effort_array_Fhistsoc <- array(NA, dim = c(length(times), length(gear_names)), dimnames = list(time = times, gear = gear_names))
effort_array_Fnat <- array(NA, dim = c(length(times), length(gear_names)), dimnames = list(time = times, gear = gear_names))

# Now fill the array
# Remember, the first 600 years are for spin-up
# During this time, the first input value is repeated at each time step
for (t in (times - 1349)) {
    if (t <= 601) { 
        effort_array_Fhistsoc[t,"Longline"] <- effort_Fhistsoc[1]
        effort_array_Fnat[t,"Longline"] <- 0
    } else {
        effort_array_Fhistsoc[t,"Longline"] <- effort_Fhistsoc[t - 600]
        effort_array_Fnat[t,"Longline"] <- 0
        }
}

Temperature

There are three climate scenarios, with different temperature forcing for each:

picontrol: Pre-industrial control
ssp126: Simulated historical climate, then SSP1-RCP2.6 climate
ssp585: Simulated historical climate, then SSP5-RCP8.5 climate

Note that only the years 1950 - 2100 are needed. Temperature forcing files were created with “PrepTemperature_GFDL.Rmd” and “PrepTemperature_IPSL.Rmd”.

# Load data for each CMIP6 model: GFDL-ESM4 and IPSL-CM6A-LR
GFDL_temperature_PIcontrol <- read.table("GFDL_ocean_temp_array_PIcontrol.dat")
GFDL_temperature_CC126 <- read.table("GFDL_ocean_temp_array_CCscenario_126.dat")
GFDL_temperature_CC585 <- read.table("GFDL_ocean_temp_array_CCscenario_585.dat")

IPSL_temperature_PIcontrol <- read.table("IPSL_ocean_temp_array_PIcontrol.dat")
IPSL_temperature_CC126 <- read.table("IPSL_ocean_temp_array_CCscenario_126.dat")
IPSL_temperature_CC585 <- read.table("IPSL_ocean_temp_array_CCscenario_585.dat")

GFDL_temperature_PIcontrol <- as(GFDL_temperature_PIcontrol, "matrix")
GFDL_temperature_CC126 <- as(GFDL_temperature_CC126, "matrix")
GFDL_temperature_CC585 <- as(GFDL_temperature_CC585, "matrix")

IPSL_temperature_PIcontrol <- as(IPSL_temperature_PIcontrol, "matrix")
IPSL_temperature_CC126 <- as(IPSL_temperature_CC126, "matrix")
IPSL_temperature_CC585 <- as(IPSL_temperature_CC585, "matrix")

# Build temperature arrays following methods used above
species <- params@species_params$species
ocean_temp_array_GFDL_PIcontrol <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))
ocean_temp_array_GFDL_CC126 <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))
ocean_temp_array_GFDL_CC585 <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))
ocean_temp_array_IPSL_PIcontrol <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))
ocean_temp_array_IPSL_CC126 <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))
ocean_temp_array_IPSL_CC585 <- array(NA, dim = c(length(times), length(species)), dimnames = list(time = times, sp = species))

for (t in (times - 1349)) {
    if (t <= 601) { 
        ocean_temp_array_GFDL_PIcontrol[t,] <- GFDL_temperature_PIcontrol[1,]
        ocean_temp_array_GFDL_CC126[t,] <- GFDL_temperature_CC126[1,]
        ocean_temp_array_GFDL_CC585[t,] <- GFDL_temperature_CC585[1,]
        ocean_temp_array_IPSL_PIcontrol[t,] <- IPSL_temperature_PIcontrol[1,]
        ocean_temp_array_IPSL_CC126[t,] <- IPSL_temperature_CC126[1,]
        ocean_temp_array_IPSL_CC585[t,] <- IPSL_temperature_CC585[1,]
    } else {
        ocean_temp_array_GFDL_PIcontrol[t,] <- GFDL_temperature_PIcontrol[t - 600,]
        ocean_temp_array_GFDL_CC126[t,] <- GFDL_temperature_CC126[t - 600,]
        ocean_temp_array_GFDL_CC585[t,] <- GFDL_temperature_CC585[t - 600,]
        ocean_temp_array_IPSL_PIcontrol[t,] <- IPSL_temperature_PIcontrol[t - 600,]
        ocean_temp_array_IPSL_CC126[t,] <- IPSL_temperature_CC126[t - 600,]
        ocean_temp_array_IPSL_CC585[t,] <- IPSL_temperature_CC585[t - 600,]
        }
}

Plankton

There is also different plankton forcing for each of the three climate scenarios. These plankton forcings are used to create the resource spectra (which include more than plankton). Plankton forcing files were created with “PrepPlankton_scaled.Rmd”.

# Load data for each CMIP6 model: GFDL-ESM4 and IPSL-CM6A-LR
GFDL_n_pp_PIcontrol <- read.table("GFDL_n_pp_array_PIcontrol_scaled_S1.0I0.85.dat")
GFDL_n_pp_CC126 <- read.table("GFDL_n_pp_array_CCscenario_126_scaled_S1.0I0.85.dat")
GFDL_n_pp_CC585 <- read.table("GFDL_n_pp_array_CCscenario_585_scaled_S1.0I0.85.dat")

IPSL_n_pp_PIcontrol <- read.table("IPSL_n_pp_array_PIcontrol_scaled_S1.0I0.85.dat")
IPSL_n_pp_CC126 <- read.table("IPSL_n_pp_array_CCscenario_126_scaled_S1.0I0.85.dat")
IPSL_n_pp_CC585 <- read.table("IPSL_n_pp_array_CCscenario_585_scaled_S1.0I0.85.dat")

GFDL_n_pp_PIcontrol <- as(GFDL_n_pp_PIcontrol, "matrix")
GFDL_n_pp_CC126 <- as(GFDL_n_pp_CC126, "matrix")
GFDL_n_pp_CC585 <- as(GFDL_n_pp_CC585, "matrix")

IPSL_n_pp_PIcontrol <- as(IPSL_n_pp_PIcontrol, "matrix")
IPSL_n_pp_CC126 <- as(IPSL_n_pp_CC126, "matrix")
IPSL_n_pp_CC585 <- as(IPSL_n_pp_CC585, "matrix")

# Build plankton arrays to be filled
sizes <- params@w_full

n_pp_array_GFDL_PIcontrol <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))
n_pp_array_GFDL_CC126 <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))
n_pp_array_GFDL_CC585 <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))
n_pp_array_IPSL_PIcontrol <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))
n_pp_array_IPSL_CC126 <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))
n_pp_array_IPSL_CC585 <- array(NA, dim = c(length(times), length(sizes)), dimnames = list(time = times, w = sizes))

# Fill arrays, remembering that values in n_pp are log10 abundances.
# They'll need to be transformed (inverse log10) and divided by bin width (dw_full)
for (t in (times - 1349)) {
    if (t <= 601) { 
        n_pp_array_GFDL_PIcontrol[t,] <- (10^(GFDL_n_pp_PIcontrol[1,]))/params@dw_full
        n_pp_array_GFDL_CC126[t,] <- (10^(GFDL_n_pp_CC126[1,]))/params@dw_full
        n_pp_array_GFDL_CC585[t,] <- (10^(GFDL_n_pp_CC585[1,]))/params@dw_full
        n_pp_array_IPSL_PIcontrol[t,] <- (10^(IPSL_n_pp_PIcontrol[1,]))/params@dw_full
        n_pp_array_IPSL_CC126[t,] <- (10^(IPSL_n_pp_CC126[1,]))/params@dw_full
        n_pp_array_IPSL_CC585[t,] <- (10^(IPSL_n_pp_CC585[1,]))/params@dw_full
    } else {
        n_pp_array_GFDL_PIcontrol[t,] <- (10^(GFDL_n_pp_PIcontrol[t - 600,]))/params@dw_full
        n_pp_array_GFDL_CC126[t,] <- (10^(GFDL_n_pp_CC126[t - 600,]))/params@dw_full
        n_pp_array_GFDL_CC585[t,] <- (10^(GFDL_n_pp_CC585[t - 600,]))/params@dw_full
        n_pp_array_IPSL_PIcontrol[t,] <- (10^(IPSL_n_pp_PIcontrol[t - 600,]))/params@dw_full
        n_pp_array_IPSL_CC126[t,] <- (10^(IPSL_n_pp_CC126[t - 600,]))/params@dw_full
        n_pp_array_IPSL_CC585[t,] <- (10^(IPSL_n_pp_CC585[t - 600,]))/params@dw_full
        }
}

Write parameters and functions for therMizer extension

The therMizer extension that allows us to use the temperature and plankton forcings includes a few new parameters that we can determine based on the user input above. There are also several functions that we’ll write, too.

Parameters

The parameter t_idx will help with the simulations by providing the correct time index for the plankton and temperature forcing during the simulations.

# Time indexing parameter
# This will be added to t to convert the year into an index for the n_pp and ocean_temp arrays
if (min(times) == 0) {
  other_params(params)$t_idx = 1
} else if (min(times) == 1) {
  other_params(params)$t_idx = 0
} else {
  other_params(params)$t_idx = -(min(times) - 1)
}

To scale the effect of temperature on encounter rate to a value ranging from 0 - 1, it is necessary to divide by the maximum possible value for each species. To scale the effect of temperature on metabolism to a value ranging from 0 - 1, it is necessary to subtract the minimum vaule for each species and then divide by the range. This requires a bit of straightforward arithmetic, and users could do this on their end if they’re so inclined. These parameters handle that math so the user doesn’t have to.

species_params(params)$encounter_scale <- rep(NA, length(params@species_params$temp_min))

for (indv in seq(1:length(params@species_params$temp_min))) {
  
  # Create a vector of all temperatures each species might encounter
  temperature <- seq(params@species_params$temp_min[indv], params@species_params$temp_max[indv], by = 0.1)
  
  # Find the maximum value of the unscaled effect of temperature on encounter rate for each species 
  species_params(params)$encounter_scale[indv] <- max((temperature) * (temperature - params@species_params$temp_min[indv]) * (params@species_params$temp_max[indv] - temperature))
}

# Determine the minimum, maximum, and range of value for the effect of temperature on metabolism
min_metab_value <- (exp(25.22 - (0.63/((8.62e-5)*(273 + params@species_params$temp_min)))))
max_metab_value <- (exp(25.22 - (0.63/((8.62e-5)*(273 + params@species_params$temp_max)))))
        
species_params(params)$metab_min <- min_metab_value
species_params(params)$metab_range <- max_metab_value - min_metab_value

Functions

Temperature will be added to the fuctions that determine encounter rate and energy available for growth and reproduction.

To scale encounter rate with temperature, we’re essentially taking a temperature-dependent proportion of the value calculated by the mizerEncounter function. If species are at their thermal optimum, we take the full value. Elsewhere in their thermal range, we take a proportion that goes to zero at the limits of species’ thermal tolerence.

therMizerEncounter <- function(params, n, n_pp, n_other, t, ...) {
  
  # Access the correct element
  temp_at_t <- params@other_params$other$ocean_temp[t + params@other_params$other$t_idx,]
  
  # Calculate unscaled temperature effect using a generic polynomial rate equation
  unscaled_temp_effect <- temp_at_t * (temp_at_t - params@species_params$temp_min) * (params@species_params$temp_max - temp_at_t)
  
  # Scale using encounter_scale parameter
  scaled_temp_effect <- unscaled_temp_effect / params@species_params$encounter_scale
  
  # Set the encounter rate to zero if temperature is outside species' thermal tolerance
  above_max <- which(temp_at_t > params@species_params$temp_max)
  below_min <- which(temp_at_t < params@species_params$temp_min)
  
  if (length(above_max) > 0)
    scaled_temp_effect[above_max] = 0
  
  if (length(below_min) > 0)
    scaled_temp_effect[below_min] = 0
  
  # Calculate maximum possible encounter rate
  max_encounter <- mizerEncounter(params, n = n, n_pp = n_pp, n_other = n_other, ...)
  
  # Apply temperature effect
  return(max_encounter * scaled_temp_effect)
      
}

To calculate the effect of temperature on metabolim, we use an Arrhenius function to scale the cost of metabolism. When species are at their thermal maximum, the cost of metabolism is at its maximum. When species are at their thermal minimum, the cost of metabolism is at its minimum

therMizerEReproAndGrowth <- function(params, n, n_pp, n_other, t, encounter,
                                 feeding_level, ...) {
  
  # Access the correct element
  temp_at_t <- params@other_params$other$ocean_temp[t + params@other_params$other$t_idx,]
  
  # Arrhenius equation
  unscaled_temp_effect <- (exp(25.22 - (0.63/((8.62e-5)*(273 + temp_at_t)))))
  
  # Arrhenius equation scaled to a value between 0 and 1
  temp_effect_metabolism <- (unscaled_temp_effect - params@species_params$metab_min) / params@species_params$metab_range
  
  # Set the EReproAndGrowth to zero if temperature is outside species' thermal tolerance
  Emultiplier <- rep(1, length(params@species_params$species))
  
  above_max <- which(temp_at_t > params@species_params$temp_max)
  below_min <- which(temp_at_t < params@species_params$temp_min)
  
  if (length(above_max) > 0)
    Emultiplier[above_max] = 0
  
  if (length(below_min) > 0)
    Emultiplier[below_min] = 0
  
  # Apply scaled Arrhenius value to metabolism
  (sweep((1 - feeding_level) * encounter, 1, params@species_params$alpha, "*", check.margin = FALSE) - params@metab*temp_effect_metabolism)*Emultiplier  
  
}

Now we’ll write a function to use the CMIP6 plankton densities rather than the mizer resource dynamics. Code for this was informed by the approach used in: https://rpubs.com/gustav/plankton-anchovy.

# Set up new resource forcing "function" - just changes to different time slot in the n_pp_array array
plankton_forcing <- function(params, t, ...) {
  return(other_params(params)$n_pp_array[t + params@other_params$other$t_idx,])  
}

Set the new rate functions and new resource

params <- setRateFunction(params, "Encounter", "therMizerEncounter")
params <- setRateFunction(params, "EReproAndGrowth", "therMizerEReproAndGrowth")

params <- setResource(params, resource_dynamics = "plankton_forcing")

Run the simulations

Because temperature and n_pp forcing are parameters, we’ll need to make unique params objects for each CMIP6 model and climate scenario. This means we’ll have six new params objects (2 models x 3 climate scenarios).

# Create parameter objects
params_GFDL_picontrol <- params
params_GFDL_ssp1rcp26 <- params
params_GFDL_ssp5rcp85 <- params

params_IPSL_picontrol <- params
params_IPSL_ssp1rcp26 <- params
params_IPSL_ssp5rcp85 <- params

# Attach temperature
other_params(params_GFDL_picontrol)$ocean_temp <- ocean_temp_array_GFDL_PIcontrol
other_params(params_GFDL_ssp1rcp26)$ocean_temp <- ocean_temp_array_GFDL_CC126
other_params(params_GFDL_ssp5rcp85)$ocean_temp <- ocean_temp_array_GFDL_CC585

other_params(params_IPSL_picontrol)$ocean_temp <- ocean_temp_array_IPSL_PIcontrol
other_params(params_IPSL_ssp1rcp26)$ocean_temp <- ocean_temp_array_IPSL_CC126
other_params(params_IPSL_ssp5rcp85)$ocean_temp <- ocean_temp_array_IPSL_CC585

# Attach plankton
other_params(params_GFDL_picontrol)$n_pp_array <- n_pp_array_GFDL_PIcontrol
other_params(params_GFDL_ssp1rcp26)$n_pp_array <- n_pp_array_GFDL_CC126
other_params(params_GFDL_ssp5rcp85)$n_pp_array <- n_pp_array_GFDL_CC585

other_params(params_IPSL_picontrol)$n_pp_array <- n_pp_array_IPSL_PIcontrol
other_params(params_IPSL_ssp1rcp26)$n_pp_array <- n_pp_array_IPSL_CC126
other_params(params_IPSL_ssp5rcp85)$n_pp_array <- n_pp_array_IPSL_CC585

Now we can run the simulations. There are 12 runs here: 3 climate scenarios x 2 fishing scenarios x 2 CMIP6 models.

When running calling project, we’ll also provide the vector for initial_n_pp, which is simply the first value from the n_pp arrays.

# GFDL-ESM4
sim_GFDL_picontrol_histsoc <- project(params_GFDL_picontrol, initial_n_pp = n_pp_array_GFDL_PIcontrol[1,], t_max = length(times), effort = effort_array_Fhistsoc)
sim_GFDL_picontrol_nat <- project(params_GFDL_picontrol, initial_n_pp = n_pp_array_GFDL_PIcontrol[1,], t_max = length(times), effort = effort_array_Fnat)

sim_GFDL_ssp1rcp26_histsoc <- project(params_GFDL_ssp1rcp26, initial_n_pp = n_pp_array_GFDL_CC126[1,], t_max = length(times), effort = effort_array_Fhistsoc)
sim_GFDL_ssp1rcp26_nat <- project(params_GFDL_ssp1rcp26, initial_n_pp = n_pp_array_GFDL_CC126[1,], t_max = length(times), effort = effort_array_Fnat)

sim_GFDL_ssp5rcp85_histsoc <- project(params_GFDL_ssp5rcp85, initial_n_pp = n_pp_array_GFDL_CC585[1,], t_max = length(times), effort = effort_array_Fhistsoc)
sim_GFDL_ssp5rcp85_nat <- project(params_GFDL_ssp5rcp85, initial_n_pp = n_pp_array_GFDL_CC585[1,], t_max = length(times), effort = effort_array_Fnat)

# IPSL-CM6A-LR
sim_IPSL_picontrol_histsoc <- project(params_IPSL_picontrol, initial_n_pp = n_pp_array_IPSL_PIcontrol[1,], t_max = length(times), effort = effort_array_Fhistsoc)
sim_IPSL_picontrol_nat <- project(params_IPSL_picontrol, initial_n_pp = n_pp_array_IPSL_PIcontrol[1,], t_max = length(times), effort = effort_array_Fnat)

sim_IPSL_ssp1rcp26_histsoc <- project(params_IPSL_ssp1rcp26, initial_n_pp = n_pp_array_IPSL_CC126[1,], t_max = length(times), effort = effort_array_Fhistsoc)
sim_IPSL_ssp1rcp26_nat <- project(params_IPSL_ssp1rcp26, initial_n_pp = n_pp_array_IPSL_CC126[1,], t_max = length(times), effort = effort_array_Fnat)

sim_IPSL_ssp5rcp85_histsoc <- project(params_IPSL_ssp5rcp85, initial_n_pp = n_pp_array_IPSL_CC585[1,], t_max = length(times), effort = effort_array_Fhistsoc)
sim_IPSL_ssp5rcp85_nat <- project(params_IPSL_ssp5rcp85, initial_n_pp = n_pp_array_IPSL_CC585[1,], t_max = length(times), effort = effort_array_Fnat)

And plot the results to get a sense of what things look like.

plot(sim_GFDL_picontrol_histsoc)

plot(sim_GFDL_picontrol_nat)

plot(sim_GFDL_ssp1rcp26_histsoc)

plot(sim_GFDL_ssp1rcp26_nat)

plot(sim_GFDL_ssp5rcp85_histsoc)

plot(sim_GFDL_ssp5rcp85_nat)

# IPSL-CM6A-LR
plot(sim_IPSL_picontrol_histsoc)

plot(sim_IPSL_picontrol_nat)

plot(sim_IPSL_ssp1rcp26_histsoc)

plot(sim_IPSL_ssp1rcp26_nat)

plot(sim_IPSL_ssp5rcp85_histsoc)

plot(sim_IPSL_ssp5rcp85_nat)

After checking through the code and results to make sure everything worked, we’ll save the sim objects so that we can prepare the output as FishMIP requests.

# GFDL-ESM4
save(sim_GFDL_picontrol_histsoc, file = "sim_GFDL_picontrol_histsoc.Rdata", ascii = TRUE)
save(sim_GFDL_picontrol_nat, file = "sim_GFDL_picontrol_nat.Rdata", ascii = TRUE)

save(sim_GFDL_ssp1rcp26_histsoc, file = "sim_GFDL_ssp1rcp26_histsoc.Rdata", ascii = TRUE)
save(sim_GFDL_ssp1rcp26_nat, file = "sim_GFDL_ssp1rcp26_nat.Rdata", ascii = TRUE)

save(sim_GFDL_ssp5rcp85_histsoc, file = "sim_GFDL_ssp5rcp85_histsoc.Rdata", ascii = TRUE)
save(sim_GFDL_ssp5rcp85_nat, file = "sim_GFDL_ssp5rcp85_nat.Rdata", ascii = TRUE)

# IPSL-CM6A-LR
save(sim_IPSL_picontrol_histsoc, file = "sim_IPSL_picontrol_histsoc.Rdata", ascii = TRUE)
save(sim_IPSL_picontrol_nat, file = "sim_IPSL_picontrol_nat.Rdata", ascii = TRUE)

save(sim_IPSL_ssp1rcp26_histsoc, file = "sim_IPSL_ssp1rcp26_histsoc.Rdata", ascii = TRUE)
save(sim_IPSL_ssp1rcp26_nat, file = "sim_IPSL_ssp1rcp26_nat.Rdata", ascii = TRUE)

save(sim_IPSL_ssp5rcp85_histsoc, file = "sim_IPSL_ssp5rcp85_histsoc.Rdata", ascii = TRUE)
save(sim_IPSL_ssp5rcp85_nat, file = "sim_IPSL_ssp5rcp85_nat.Rdata", ascii = TRUE)

run_therMizer_FishMIP

Phoebe.Woodworth-Jefcoats@noaa.gov

8/25/2020